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ABSTRACT 

The  purpose  of  this  report  is  to  discuss  one  parameter 
methods  for  the  numerical  integration  of  stiff  differential  equations, 
i.e.,  systems  with  widely  separated  time  constants.   The  methods  which 
are  optimal  in  terms  of  least  amount  of  solution  history  data  to  "be 
saved  at  each  step  of  the  corrector  iteration  have  been  known  for  order 
as  high  as  six.   We  find  that  in  one  parameter  methods  the  stability 
and  accuracy  are  controlled  by  the  same  variable  y.      Suitable  values 
of  y  lead  to  methods  which  have  either  optimal  accuracy  or  optimal 
stability.   The  methods  of  order  six  have  been  presented.   The  section 
of  the  locus  p(£)/o(0  for  methods  of  order  six  is  shown  in  Figure  7. 
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6.  SIXTH  ORDER  FORMULA 16 
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1.      INTRODUCTION 

Many  numerical  problems   in   applied  mathematics   entail  the 
solution  of  systems   of  ordinary   differential  equations  with  a  property 
given  by  the  following  definition. 

Definition.     A  system  of  ordinary  differential  equations  y'  =  f(ty   y)3 
y(a)  =  n  is  said  to  be  stiff  if  the  eigen  values  of  the  matrix  (-^-)  has 
at  every  -point  t  negative  real  parts,   and  differ  greatly  in  magnitude. 

Equations   of  this    form  frequently   arise   in  physical  equilibrium 
problems   such   as   chemical  or  nuclear  reactions  where   there   are  many 
components  with  a  large   range   of  time   constants.      The   standard  numerical 
integration  methods   are   unstable   if  the   step  size   used  is  much   greater 
than  the   smallest   time   constant.      Dahlquist    [l]   has   introduced  the 
concept   of  A-stability.      The   theorem  of  Dahlquist   states   that   among  all 
linear  multistep  methods   the   trapezoidal  rule   is   the  most   accurate. 
Gear    [2,3]  has    developed  stiffly   stable  methods  based  on  the  necessary 
requirements   of  the   stiff   differential  equations.      He  has   obtained  stiffly 
stable  methods   of  order  as   high   as   six.      Dill    [k]   has   extended  this 
analysis   to  the  methods   of  order  seven   and  eight.      Jain   and  Srivastava   [5] 
examined  a  class   of  methods   and  obtained  stiffly  stable  methods   of  order 
eleven.      The  purpose   of  this   report   is   to  obtain   optimal  stiffly  stable 
methods   for  the   stiff   differential  equation 


(1) 


|r=  f(t,  y),  y(a)  =  n 


2.   OPTIMAL  STIFFLY  STABLE  METHODS 

The  requirements  for  the  stiffly  stable  methods  are  shown 
in  Figure  1. 


STABLE 


ACCURATE 


S   T  A  B   L  E 


•0 


hX  PLANE 


Fi  gure  1 
STABILITY  AND  ACCURACY  REGIONS 
The  existence  of  stiffly  stable  methods  depends  on  the 
parameters  D,  0,  a  and  on  the  definition  of  accuracy.   The  usual  definition 
of  accuracy  is  that  of  the  order.   Dahlquist  has  shown  that  if  D  is  zero 
then  the  order  of  the  method  cannot  exceed  two.   Here  we  have  obtained 
stiffly  stable  methods  with  extended  region  of  stability  of  order  as  high 
as  six  for  suitable  parameters  D,  0,  and  a. 


Third  Order  Formula:   Multistep  formulas  may  be  obtained  by  writing  a 
linear  relation  between  values  of  the  function  and  derivative  at  specified 
points  and  determining  the  coefficients  in  the  relation  by  expansions  in 
Taylor  series  about  an  arbitrary  point.   For  instance,  we  may  write 

(2)  y  j.n  =  a-,  y    +  ao  y    -,  +  a->  y    «  +  ai  y    o  +  h$o  y'^n  +  T 

n+l  In         2  "n-1  3     n-2  k     n-3  0     n+l         n 

The   fourth  order  formula,   that  is,   the  one  whose  truncation  error  is   of 
order  h    ,   has  been   obtained  by  Gear 

(3)  =  U8     36     +16       _I     +  h  12  . 

yn+l  "  25  yn   25  yn-l   25  yn-2   25  yn-3     25  yn+l 

12  ,5   (5) 
"125h  yn 

If  instead  of  requiring  that   T    ,   the  truncation  error,    should  be   of 

fifth  order,   we   specify  a  fourth   order  truncation  error  we   find  that 

there   is    a  one-parameter  family  of  methods. 

This   may  be  written 

(h)  .   ,U8  _  _26_  )  , 36  _  _5T    v  ,16        ^2    > 

yn+l  "    V25   "   300Y;yn        ^25        300Y;   yn-l        ^25   "   300Y;   yn-2 

-    (—  -  —  v)y  +  h(—  +  —  y)   v'         +   T 

v25        300r/,yn-3  25        300 w   °n+l         n 

or  alternatively 

1+8  36  16  3  .12     . 

v         =  —  v     -  —  v         +  —  v         -  —  v         +  h  —  v 
°n+l       25  "n       25  °n-l       25  °n-2       25  ^n-3  25     n+l 

[26y^    -   5Ty^    n    +    ^    _    -   lly      _   -  h6y «    .  ]   +  T 


where 


300  'n        "Jn-1  °n-2  Jn-3  "'n+l  n 


YhVU)/: 

m     _  i  n 


{ 


n  5    (S) 

c   12b/y^V125 


This  later  way  of  writing  emphasizes  the  nature  of  the  integration 

formula.   It  is  seen  to  be  composed  of  the  fourth  order  formula  with 

the  minus  of  a  third  order  restrictive  condition.   Note  also  that 

y  =  —   gives  the  third  order  stiffly  stable  method 

18      9        2      L  ,  6   ,      3  A   (U) 
yn+l  =  iT  yn  "  IT  yn-l  +  IT  yn-2  +  h  IT  yn+l  "  22  h  yn 

The  discretization  error  of  the  formula  {h)    is  defined  as  the  difference 

between  the  value  y   calculated  from  (k)    and  the  exact  solution  y(x  ). 

n  n 

Define  e      =  y      -  y(x   ) 

n  n  n 

Then  the   error  e      obeys   the    difference   equation 

(5)  (iii      _2£    \  (36       _5T    n  / 16         k2     > 
n+1  "   ^25   "   300Y;£n   "   ^25   "   300YJ£n-l  +   l25   "   300Yj£n-2 

-(2T-3^Y)£n-3  +  hA(i+34Y)£n+l-Tn 
for  the    differential  equation  y'    =   Ay. 

This   linear  difference  equation  with   constant   coefficients  may  be  solved 
by  setting  e     =   e,    ,  which  gives  the   characteristic  equation 

(6)  pU)    -  hA  cr(5)   +  T     =0 

n 


where 


and 


PU;        *  l25        300YH     +   l25        300YJ^     "    l25        300YH 

v25        300 w 


ff(€)    =    (2?+    300Yk 


It  is  necessary  to  bound  the  solution  of  the  inhomogeneous  difference 

equation  (5),  Henrici  [6],   It  depends  on  the  stability  of  the 

corresponding  homogeneous  equation  which  is  obtained  when  T  is  set 

to  zero.   The  difference  equation  is  stable  if  and  only  if  all  roots  of 

the  polynomial  equation 

(T)        pU)  -  hX  c(c)  *   0 

are  inside  the  unit  circle  or  on  the  unit  circle  and  simple.   If  stiff 

differential  equations  are  to  be  integrated  with  large  values  of  h,  then 

large  values  of  hX  must  not  make  (7)  unstable.   Letting  hX  ->  °°,  the 

roots  of  (7)  approach  those  of  a(C)  =  0.   This  implies  that  the  polynomial 

a(£)  must  not  have  roots  outside  the  unit  circle  and  those  roots  £.  for 

l 

which    |£.|    =  1  are   simple.      This    condition   is   already  satisfied  by  a(£) 
in  this    case.      For  stiff  stability  we  want  hX  values   such  that    (7)   has 
roots   inside  the   unit    circle   or  on  the  unit    circle   and  simple.      This 
region  is  bounded  by  the   locus   of  p(0/a(0    in  the  hX-plane   for   £  =  e      , 
0e[O,   2tt].      For  different   values   of  y  the   locus   of  p(£)/a(£)    is   shown 
in  Figure   2. 

The   locus  boundaries   are   closed  and  symmetrical  about  the   real  axis, 
pass   through  the   origin   and  steadily   decrease   in   area  as  y   increases   from 
0  to   7,    approximately.      At  y   =   7»   the   formula  ceases   to  be   stiffly  stable. 
The   intercept  of  the  boundary  on  the   real  axis   is    given  by 

(8)  ™   -  k  192_-_lJx 

hA    "    3      2U   +      y 

36  . 
y  =  —  gives  third  order  formula  which  is  optimal  in  terms  of  least 

amount  of  solution  history  data  to  be  saved  at  each  step  of  the 

corrector  iteration. 


Here  we  find  that  the  values  of  y  in  the  interval 

0  <  Y   <  ■=—  give  third  order  formulas  which  are  stiffly  stable  with 

reduced  truncation  error  coefficient.  We  shall  denote  such  formulas  by 

P.   Similarly  for  the  values  of  Y  in  the  interval  —  <  y  <  7,  we  obtain 

third  order  stiffly  stable  formulas  which  have  extended  region  of 

stability.   We  shall  denote  such  formulas  by  Q.   The  locus  of 

P(0/a(C),  £  =  e1  ,  0e[O,  2tt  ]  for  y  =  6  is  shown  in  Figure  3.   The  P 

and  Q,  type  formulas  for  y  =  2  and  y  =  6  are  listed  here.   The  values  of 

parameters  C  (error  constant),  D,  0  and  max|?|  (£.  being  the  roots  of 

pCO  =  0)  for  formulas  of  order  p  ^_  6  are  listed  in  the  table. 

(P)  .262      159       _5jt       _L_      1  h  T8   ' 

yn+l  "   150  yn  "  150  yn-l   150  yn-2  "  150  yn-3     150  yn+l 

12   Jn 

(Q)  _  1U      _3_       _2_       _1_      +hA. 

yn+l  "  10  yn  "  10  yn-l  "  10  yn-2   10  yn-3     10  yn+l 

1  h    (k) 

One  parameter  family  of  formulas  together  with  P  and  Q  type  formulas 
for  order  six  are  given  here. 


Fourth  Order  Formula:   One  parameter  family  of  fourth  order  formula 

is   given  by 

(9)  .   300_  300  200_  _J_5_  _12_ 

yn+l  ""  137  yn  "  137  yn-l       137  yn-2   "   137  yn-3       137  yn-4 

+  h  §fK+l  '   3288   [TTyn   "  21K-1  +  23K-2 
~  122yn_3  +  25yn_u  -  h  I2yi+1]   -  ^  h5y [ 5 } 


The   characteristic  equation   is 

<»>  C5  -  ^  (300  -  '%,)+  ♦  ^  (300  -  $,H3  -  ^  (200  -  f  YH2 

+  I5r  (75  -  W^  -  I5T  (12  -  1^  "  hx  iff  (60  +  W^  "  ° 

The  locus  of  p(£)/a(0  for  different  values  of  y  is  shown  in  Figure  k. 

pQO 

Y  =  ~  gives  fourth  order  formula  obtained  by  Gear,  y   >_  25   does  not 

give  stiffly  stable  methods.   The  formulas  for  y  =   5  and  y   =  2k   are 

listed  below. 

(P)  6815      6130        3630        1190        163 

yn+l  "  3288  yn  "  3288  yn-l  +  3288  yn-2  "  3288  yn-3   3288  Yn-k 

l  u  J^OQ  1     x  1,5  (5) 
+  h  3288  yn+l  "  2¥  h  yn 

(Q)  .  223     _86_      J}L  +  J±l.  13 

yn+l  ""  137  yn   137  yn-l  "  137  yn-2   137  yn-3   137  Vn-k 


Fifth  Order  Formula:   One  parameter  family  of  fifth  order  formula  is 

(H)  1        l^r,  522     v  1        ,,„  1753     N 

yn+l  =  IW  (36°   "   725Y)yn  "  HPf  {k5°  ~  l2TY)yn-l 

1        n  no  25U0      N  1         ,ooc  1980     s 

+   1U7   (U0°   "  ^^Y)yn-2   "  IW  (225   -  72TY)yn-3 


1     ,VQ       810    .  _1_  (  137   x 

1U7  "   720^11-1*  "  lU7   U     "   720Y;yn-5 


6  (6) 


+  h 


lU7  720  wyn+l  720 


The   locus   of  p(£)/a(£)    for  different  values   of  y   is   given   in  Figure   5. 

7200      . 
Y  -     -,  o7    gives    fifth  order  stiffly  stable   formula  and  has   the  virtue 

of  economy  of  storage.      For  y   >   HO,    fifth   order  stiffly  stable    formula 

ceases   to  exist.      The   stiffly  stable   formulas   for  y  =   36   and  y  =   96  becomes 


(p)          .  6678  72U5  5^6o      2520 

yn+l  "  29^0  yn  29^0  yn-l  29^0  yn-2   29^0  yn-3 

630  63  ,    1260   ,     _i_h6(6) 

29U0  yn-l+  "  29^0  yn-5  29^0  yn+l  "  20   yn 

(Q)          ,  U356  32  Uo  ,  920     +  _585_ 

yn+l  "  2205  n  ~  2205  yn-l  2205  Yn-2   2205  Yn-3 

5hO  22k  L  .  1020   ,      2   6  (6) 


2205  Jn-U   2205  ,xn-5     2205  Jn+1   15   Jn 


Sixth  Order  Formula: 


(12)  y         =  -1_  (29U0   -  ^Y)V     -  -A-  (MlO  -  ^Y )y 

yn+l       1089    iy  720YJyn       1089^  720  ' yyn-l 

•       1      (U9oo-il£yW         -      1      (3675  -^Y)y 


1089  720T/Jn-2        1089  720  "'n-3 

+  IbV  (1T6U  -  W*)yn-U  "  1559    ik9°  ~  W^n-5 

+  im  (6°   -  #)yn-6  +  h  IW  (U2°  +  W)yn+1 

Y      h7  (7) 
"  3W  h  yn 


U3200 


Y  =      1  k 7     gives   six  order   stiffly   stable    formula  as   obtained  by   Gear. 

Y  =   0   gives   a  seventh   order  formula  which  is   unstable.      The   locus   of 

p(£)/a(0    as   a  function   of  y   is   shown  in  Figure   6.      The   stiffly  stable 

formulas    for  y  =  2^0   and  y  =   360   are   listed  here.      The   formula  is 

stiffly  stable   in  the   interval  15    <  Y   K  620,    approximately. 

(P)  _   8151  10593  9955  6105 

yn+l  "   3267  yn  "    3267    yn-l  '    326>7~  yn-2   "   32oT  7n-3 

2277  U51  33  1320      , 

3267  yn-U  "   3267  yn-5        3267  yn-6    "        3267  yn+l 

21       Jn 

(Q)  .   1737  2061  1685  810 

yn+l  "      726    yn         72  6   yn-l         726    yn-2   ~   726  yn-3 
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1   ,7   (7) 
"  IThyn 


TABLE  OF  VALUES  OF  THE  PARAMETERS 


Order  and  Type 

of  the  Method    Y 


max I  £ | 
5*1 


3 

P 

■  '  ■ 
2 

-0.1777 

8U°lU' 

O.U3635 

-25/156 

Q 

'     6 

-0.0535 

86°50' 

O.I49 

-5/12 

it 

P 

5 

-1.3953 

69°20' 

0.62757 

-137/1500 

Q 

2h 

-0.U238 

77°26' 

0.60328 

-137/360 

5 

P 

36 

-3.2166 

59°12' 

0.7^3 

-1U7/1260 

Q 

96 

-I.U966 

66°3' 

0.71381 

-1U7/510 

6 

p 

2U0 

-7.1812 

U8°36» 

0.88UUU 

-363/3080 

Q 

360 

-5.OIU9 

50°3)+' 

0.8^820 

-121/700 

10 

3.   CONCLUSIONS 

The  stiffly  stable  methods  have  been  found  very  efficient 
in  integrating  stiff  differential  equations  because  they  allow  a  much 
larger  step  size  while  maintaining  stability.   Here  we  have  developed 
one  parameter  family  of  methods  suited  to  integrating  stiff  differential 
equations.   It  is  found  that  stability  and  accuracy  are  controlled  by 
the  same  variable  y   Ihe  methods  which  are  optimal  in  terms  of  least 
amount  of  solution  history  data  to  be  saved  at  each  step  were  obtained 
by  Gear  for  order  as  high  as  six.   A  suitable  value  of  y  leads  to 
methods  which  have  optimal  accuracy  or  optimal  stability  without 
violating  the  condition  of  economy  of  storages.   There  is  a  great 
advantage  to  be  gained  by  the  use  of  a  formula  with  optimal  stability 
when  solving  sets  of  simultaneous  differential  equations  whose  solution 
approximate  to  exponentials  with  real  negative  arguments.   An  optimal 
stability  enables  the  interval  of  integration  to  be  increased  without 
causing  instability.   However,  in  general  the  value  of  y  will  vary  with 
the  behavior  of  the  equations  we  are  attempting  to  solve.   Parametric 
methods  for  the  solution  of  ordinary  differential  equations  have  also  been 
discussed  by  Hamming  [7]  and  Robertson  [8].   Applications  of  these  and 
similar  formulas  to  predictor-corrector  methods  will  be  discussed  separately, 
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Figure  2 
Third  Order  Formula 
STIFFLY  STABLE  REGION  AS  A  FUNCTION  OF  y 
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Figure    3 
Third  Order  Formula 
LOCUS   OF  p(5)/a(5),    i  =   e10,    0e[0,    2tt] 
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Figure  h 
Fourth  Order  Formula 
STIFFLY  STABLE  REGION  AS  A  FUNCTION  OF  y 
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Figure    5 
Fifth   Order  Formula 
STIFFLY    STABLE   REGION   AS   A  FUNCTION   OF  y 
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Figure  6 
Sixth  Order  Formula 
STIFFLY  STABLE  REGION  AS  A  FUNCTION  OF  y 
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SECTION  OF  LOCUS  OF  pU)/o(£)  FOR  p  <  6 
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